%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This Matlab program is for 'The explicit formula for the Hodrick-Prescott filter in finite sample' by Adriana Cornea-Madeira
%This probram gives an alternative Matlab code to the code from HPweightsvsHPinvmatrix.m corresponding to Corollaries 
%6 and 7 which works faster than HPweightsvsHPinvmatrix.m for large sample size (denoted n). The code is for n an even number.
clear all; clc;

alpha = 1600; % smoothness parameter
l=2; % l = 2 because n is even (see (270)
n = 1000; % sample size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HP filter weights computed using the formulae from Corollaries 6 and 7, 
% for sample size even and odd. This code is faster for large n.

tic
% n even
m = n/2;
fun = @(x,z) x.*z;
i = 1:m;
x = (-1).^(i-1);
z = ones(m,1);
J = bsxfun(fun,x,z);

fun = @(i,j) sqrt(2/(n+1))*sin(i.*j);
j = i'*pi/(n+1);
T1 = bsxfun(fun,i,j);

fun = @(i,j) sqrt(2/(n+1))*sin(i.*j*2);
M2 = bsxfun(fun,i,j);

i = i';
j = i;

a1 = sin((2*j)*pi/(n+1));
b1 = sin(2*(2*j)*pi/(n+1));
c1 = 1+alpha*(2-2*cos((2*j)*pi/(n+1))).^2;

fun = @(a,b) (2/(n+1))*a.*b;
a = (2*a1-b1)./c1;
b = a';
G2 = bsxfun(fun,a,b);

k2 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*a1-b1).^2./c1));
E = k2*M2*G2*M2';

% eigenvalues
b = cos(j*pi/(n+1));
gam1 = 1./(1+4*alpha*(1-b).^2);
gam2 = 1./(1+4*alpha*(1+flipud(b)).^2);

X = T1.*J; 

V1 = T1*bsxfun(@times,gam1,T1) + X'*bsxfun(@times,flipud(gam2),X);
V2 = T1*bsxfun(@times,gam1,fliplr(X'))+(-1)^l*X'*bsxfun(@times,flipud(gam2),fliplr(X')).*J;
j = i-1;
a2 = sin((2*j+1)*pi/(n+1));
b2 = sin(2*(2*j+1)*pi/(n+1));
c2 = 1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2;
k1 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*a2-b2).^2./c2));
    
fun = @(i,jj) sqrt(2/(n+1))*sin(i.*jj);
jj = (2*j+1)'*pi/(n+1);
M1 = bsxfun(fun,i,jj);
fun = @(a,b) (2/(n+1))*a.*b;
a = (2*a2-b2)./c2;
b = a';
G1 = bsxfun(fun,a,b);
  
D = k1*M1*G1*M1';
     
H = V1+D+E;
DE = D-E;
    
% the inverse
Finv_new = [H V2+fliplr(DE);rot90(V2,2)+flipud(DE') rot90(H,2)]; 
toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare the computation time of the above code with the following code
% which also corresponds to Corrolaries 6 and 7
tic
n = 1000;
m = floor(n/2); 
% permutation matrix
p = m:-1:1;
Pm = eye(m);
Pm = Pm(p,:);
    
i = (1:1:m)';
j = i;
T1 = sqrt(2/(n+1))*sin(i*j'*pi/(n+1)); % first block from the matrix of eigenvectors
x = repmat(i',m,1);
J = ones(m,m);
J = J.*(-1).^(x-1); 
G2_0 = (2*sin((2*j)*pi/(n+1))-sin(2*(2*j)*pi/(n+1)));
Ei = (2-2*cos((2*j)*pi/(n+1))).^2;
M2 = sqrt(2/(n+1))*sin(i*(2*j)'*pi/(n+1)); 
b = cos(j*pi/(n+1));
T11 = (T1.*J)'*Pm;
T12 = Pm*(T1.*J)'*(Pm.*J);
T13 = Pm*T1.*J;
    
j = i-1;
G1_0 = (2*sin((2*j+1)*pi/(n+1))-sin(2*(2*j+1)*pi/(n+1)));
Eii = (2-2*cos((2*j+1)*pi/(n+1))).^2;
M1 = sqrt(2/(n+1))*sin(i*(2*j+1)'*pi/(n+1)); 

G2 = (2/(n+1))*(G2_0./(1+alpha*Ei))*(G2_0./(1+alpha*Ei))';
k2 = 2*alpha/(1-2*alpha*sum((2/(n+1))*G2_0.^2./(1+alpha*Ei)));
E = k2*M2*G2*M2'; 
    
% eigenvalues
b1 = (1-b);
gam1 = 1./(1+4*alpha*b1.*b1);
Pmb = (1+Pm*b);
gam2 = 1./(1+4*alpha*Pmb.*Pmb); %if n odd, the (m+1)'s eigenval is skipped but we do not need it here.

Dgam1 = diag(gam1); Dgam2 = diag(gam2); 
V1 = T1*Dgam1*T1 + T11*Dgam2*T13;
V2 = T1*Dgam1*T11 - T11*Dgam2*T12;
    
% n is even
k1 = 2*alpha/(1-2*alpha*sum((2/(n+1))*G1_0.^2./(1+alpha*Eii)));
G1 = (2/(n+1))*(G1_0./(1+alpha*Eii))*(G1_0./(1+alpha*Eii))';
D = k1*M1*G1*M1';
H = V1+D+E;
   
% the inverse
Finv = [H V2+(D-E)*Pm;Pm*(V2*Pm+D'-E') Pm*H*Pm]; 
toc